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ABSTRACT 


A recombination mechanism that significantly reduces charge 
collection from very dense ion tracks in silicon devices was 
postulated by Zoutendyk et al. The theoretical analysis presented 
here concludes that Auger recombination is such a mechanism and 
is of marginal importance for tracks produced by 270-MeV krypton, 
but of major importance for higher density tracks. The analysis 
shows that recombination loss is profoundly affected by track 
diffusion. As the track diffuses, the density and recombination 
rate decrease so fast that the linear density (number of elec- 
tron-hole pairs per unit length) approaches a non-zero limiting 
value as t-»<». Furthermore, the linear density is very nearly 
equal to this limiting value in a few picoseconds or less. When 
Auger recombination accompanies charge transport processes that 
have much longer time scales, it can be simulated by assigning a 
reduced linear energy transfer to the ion. 
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1. INTRODUCTION 


Zoutendyk et al. [1] postulated the existence of a recombina- 
tion mechanism that significantly reduces charge collection from 
very dense ion tracks in silicon devices. This can impact device 
susceptibility to single-event effects (SEE) . It is reasonable to 
ask if Auger recombination (AR) can have such an effect. This 
report presents a theoretical investigation of the impact of AR 
on charge collection from an ion track. 


2. PHYSICAL MODELS 


Ion tracks consist of high electron-hole pair (EHP) density, 
which is well described by the ambipolar diffusion equation 

D div grad P - SP/St = R (1) 

where D is the ambipolar diffusion constant and P is the excess 
hole density (hole density minus equilibrium hole density) . 
Because of quasi-neutrality and high carrier density, P is also 
equal to the EHP density. R is the recombination rate and S/St 
denotes partial derivative with respect to time. It will be seen 
that the time scale for significant AR is a few picoseconds or 
less, which is very short compared to typical lifetimes for 
Shockley-Read-Hall recombination (SRHR) . Therefore, SRHR has a 
negligible effect on EHP loss via AR, and the objective is to 
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estimate only this loss, so only AR is considered in this analy- 
sis. For quasi-neutral and high-density conditions, the AR rate 
used by the numerical simulation code PISCES [2] reduces to 

R = CP 3 * * 


where C=C n +C p with C n =2.8xl0 -31 cm 6 /s and C p =9.9xl0~ 32 cm 6 /s 
(default values used by PISCES) . Therefore, the governing equa- 
tions are 


D div grad P - $P/6t = CP 3 for t > 0 (2) 

P = Pj for t = 0 ( 3 ) 

where Pj is the initial EHP density. The track is treated as 
infinitely long and uniform in the vertical (z) coordinate so 
that P ( r , t ) depends only on the radial distance r (where 
r 2 =x 2 +y 2 ) and on time t. The objective is to obtain an estimate 
of P and of the linear density (number of EHPs per unit length) , 
which is the x and y integral of P. 


3. MATHEMATICAL ANALYSIS 

The first objective is to solve (2) and (3) for P. These equa- 

tions will be converted into an integral equation, but the proofs 

of some theorems to follow require the integral equation to be 
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solvable by iteration (i.e., iteration produces a converging 
sequence of iterates) . To accomplish this, (2) will be rewritten 
so that the right side has a bounded P derivative. Define P M by 

P M = max Pj(r). (4) 

It is evident that P(r,t)<P M for all r,t and, therefore, the 
solution P is not changed if (2) is replaced with 

D div grad P - <5P/6t = CF(P) (5) 

where F is defined by 


r 3 if I r | < P M 

P M 2 r if | r | > P M . 


( 6 ) 


Although the right side of (5) is indistinguishable from the 
right side of (2) when P is the correct solution, there can be a 
difference in the equations when iteration is used to solve them 
because some initial guesses can cause the equations to produce 
different sequences of iterates. The motivation for (5) is to 
make convergence proofs easier. 


(A) Integral Equations 


Define the constant B by 


B ■ 3CP m 2 


(7) 
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and write (5) as 


div grad P - D -1 (BP + <5P/6t) = D -1 [CF(P) - BP]. (8) 

Define Pq as the solution to 

div grad P Q - D _1 <SP 0 /<St = 0 for t > 0 (9) 

p 0 = p i for t = 0 (10) 

and define P B by 

P B s e_Bt p 0 ( 11 ) 

so that P B satisfies 

div grad P B - D -1 (BP B + SP^St) = 0 for t > 0 (12) 

P B = P I for t = 0. (13) 


The two-dimensional infinite-space Green's function for the 
operator on the left of (9) is [3] 

G(X,t;X',t') = (t-t 1 ) -1 exp[-{ (x-x') 2 +(y-y') 2 )/{4D(t-t') }] (14) 
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where uppercase X denotes the vector (or point in space) with 
components x and y. G can be used to solve (9) and (10) for P 0 
with the result [3] 

% 

p 0 = (47TD)" 1 P I (r») G(X,t;X',0) d 2 X' (15) 

where the integral is defined by 

% P 00 f 00 

d 2 X' = dx' dy ' . (16) 

—00 —00 

J si J 

The B on the left of (8) can be eliminated by a change in vari- 
ables produced by multiplying P by e B *" and the equation can then 
be inverted by using G [3]. The final result is 

't r 

P = P + (4ttD) -1 [BP-CF (P) ] e _B(t_t,) G d 2 X’ dt ' . (17) 

JO J 

v 

For notational brevity, the arguments of P and G in the integrand 
are not displayed, but P in the integrand is P(r',t') and G is 
G (X, t ;X 1 , t ' ) . The analysis leading to (17) also applies if B is 
replaced with zero and the result is 

rt r 

p = p - C ( 47rD) -1 F (P) G d 2 X 1 dt'. (18) 

JO J 
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(B) Bounds by Iteration 


The method used in later sections to estimate P and the linear 
density will consist of constructing lower and upper bounds. Some 
concepts, discussed in this section, will be needed for future 
use. 

We will call a function P-^ an upper bound for another function 


P-j^rjt) > P 2 (r,t) for all r,t > 0. (19) 

P x is called a lower bound for P 2 if and only if P 2 is an upper 
bound for P^ 

Let p( 2 ) be obtained by starting with some pC 1 ) (parentheses 
emphasize that the superscripts are not exponents) and iterating 
(17); i.e., they are related by 

[B p( 1 ) -CFfP^ 1 ) ) ] e ”B(t-t 1 d 2 X' dt'. (20) 

It is easy to show that Bt-CF(t) is an increasing function of t 
and, therefore, if P^ 1 ) is an upper bound for P, the right side 
of (20) is an upper bound for the right side of (17), i.e., p( 2 ) 
(the left side of (20)) is an upper bound for P (the left side of 
(17)). Therefore, if we iterate (17) by starting with an upper 
bound for P, another upper bound is obtained. Repeated iteration 


P^ 2 ) = p B + ( 47rD ) -1 
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produces a sequence of upper bounds. Similarly, lower bounds 
produce new lower bounds. Equation (18) has a similar property 
except that the negative sign in front of the integral causes 
upper bounds to produce lower bounds and vice versa. It is shown 
in Appendix A that a sequence of functions obtained by iterating 
(17) converges, and if and p( 2 ) are related by (20) , then 
p(2) will be a better estimate of P than p^ 1 ) is (unless p( 1 )=p, 
in which case p( 2 )=P) 1 . Equation (18) does not have this property 
and iteration typically produces a wild oscillation (in the sense 
that if the initial function is a little larger than P, the first 
iterate will be much smaller than P and the second iterate will 
be much much larger than P, etc.). However, (18) will produce a 
good (accurate) bound if we start with a sufficiently good bound. 

The convergence proof (Appendix A) applies only for finite t. 
Functions produced by iteration of (17) have an unfortunate 
property in the limit as t^°o. This can be seen by integrating 
(20) while using 

% 

G d 2 X = 4 ttD ( 21 ) 


to get 


1. As an incidental point, if B was 3CP M /2 instead of 3CP« , 
the iteration would converge faster. But (18) would lose the 
property of bounds producing new bounds. Without this property, 
the only simple error estimates available are those derivable 
from Banach's fixed point theorem, which are extremely pessimis 

tic. 
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dt' . 


P 2 d 2 X 


P B d 2 X + 


't 

0 

V 


(BPj^-CP;,^ 3 ) d 2 X' 


e ~B (t-t 1 ) 


For large t, the first integral on the right goes to zero and, 
because of the exponential coefficient, only the asymptotic (in 
t') value of 


(BPj-CPjl 3 ) d 2 X' 


contributes to the t' integral. Assuming that the P-^ integral 
goes to zero, we have 


lim 


% 


■> 


p 2 d 2 x = 

lim 

(BP 1 ~CP 1 3 ) d 2 X' 


V 

t 1 -«» 

% 


% 


rt 


e-B(t-t') dt , 


lim 

t-*<» 


P n d z X. 


Therefore, if we define the linear densities N-^t) and N 2 (t) by 


Ni(t) = 


Pj_ d 2 X for i = 1,2 


we conclude that N 2 has the same limiting value as N-j^. In a 
sequence of iterates obtained by iterating (17), each function 
produces the same limiting value for the linear density as the 
initial guess. Fortunately, this problem is not too serious when 
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working with upper bounds for P. The reason is that if we start 
with an upper bound p(°) and then iterate (17) n times to obtain 
an upper bound P< n > , we will find that the linear density N< n ) 
first decreases and then increases in t. But the actual linear 
density is strictly decreasing so an upper bound at a given time 
is also an upper bound for all later times. Therefore, must 
be plotted in t only until it begins to increase. The minimum 
value can be used for all later times. 

The above statements indicate that the following methods will 
be effective for constructing bounds for the linear density: 

(i) To obtain an upper bound: 

1: Start with an upper bound p(°) for P and iterate 
(17) n times to obtain an upper bound p( n ) with corre- 
sponding linear density N< n >(t). Plot N< n > against t 
until it reaches a minimum and use the minimum value 
for all later times 

or 

2: Start with a lower bound f or p and iterate (17) 
n times to obtain a lower bound P^ n ^ . Then iterate (18) 
once to obtain an upper bound p( n+1 ) . 

(ii) To obtain a lower bound: 

1: Start with a lower bound P^) for P and be content 
with it (iterating (17) will not improve the limiting 
value of the linear density) 
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or 


2: Start with an upper bound p(°) for P and iterate 
(17) n times to obtain an upper bound p( n ) . Then iter- 
ate (18) once to obtain a lower bound p( n+ l) . 


(C) An Initial Estimate 

P can be estimated by constructing upper and lower bounds and 
the previous discussion shows how bounds can be constructed. But 
those methods require that we start with a bound and another 
method, discussed here, is needed to construct an initial bound. 
The idea is to guess at a function and then perform a test to 
determine if the function is a bound for P. One such test is the 
following. Let P £ be a bounded and well behaved (in the sense 
that div grad P E and $P E /<5t exist) function that is to be tested. 
P E is required to satisfy 


P E = Pj at t = 0. 


( 22 ) 


Define g by 


g = -D div grad P £ + <SP E /<St + CP E 3 . (23) 

It is shown in Appendix B that: 
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If g (r , t) >0 for all r,t>0, then P E is an upper bound for P. (24) 


If g(r, t) <0 for all r,t>0, then P £ is a lower bound for P. (25) 

If neither (24) or (25) are satisfied, the test is inconclusive 
and we should guess at another P E and try the test again. 

(D) A Particular Estimate when P j is a Gaussian Function 

Consider the special case 

P I (r) = N 0 (47rDt 0 ) _1 exp[-(4Dt 0 ) _1 r 2 ] (26) 

for some constants Nq (the initial linear density) and tg. The 
form (26) is not a good approximation for an actual ion track but 
it greatly simplifies the analysis and can produce useful esti- 
mates when combined with judicious curve fitting (this is illus- 
trated in the next section) . The integral in (15) can be evaluat- 
ed with the result 

P Q (r, t) = N 0 [47TD(t+t 0 ) ] -1 exp{-[4D(t+t 0 ) ] _1 r 2 ) . (27) 

Define p ei b y 


P E1 ^ [ 2Ct + (f(t)P 0 ) _2 ] _1 / 2 


(28) 


where f satisfies 
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f (0) = 1 


(29) 


f(t) > 0 for all t > 0 (30) 
but otherwise remains to be determined. 

Using the chain rule to evaluate the derivatives on the right 
side of (23) , we find that several terms subtract out (which is 
the motivation behind the particular estimate P E1 ) , and the 
result is 


g = fP 0 (2Ctf 2 P 0 2 + 1)-V2 gi (31) 

where g^ is defined by 

g x a (2Ctf 2 P 0 2 +l)f '/f +[3C/(2D) ]t(t+t 0 ) -2 f 2 P 0 2 r 2 (32) 

and the prime denotes differentiation. P E1 will be constructed to 
be a lower bound for P by selecting f to make g satisfy (25) . 
This can be accomplished by selecting f to make g^^ satisfy 

max g x (r, t) = 0 for all t > 0. (33) 

r>0 

Using (27) , we find that the maximum in r is a relative maximum 
given by 
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max g^iTft) = 
r>0 

(f '/f)+3C(47rD)" 2 t(t+t 0 )“ 3 N 0 2 f 2 exp[(2/3) (t+t Q ) (f'/f)“l] 


so (33) becomes 


f'/f = — ( 3/ 2 ) Kt (t+t Q ) 3 f 2 exp[ (2/3) (t+t Q ) f'/f ] (34) 


where 


K s 2CN Q 2 e 1 ( 4 ttD) 2 . (35) 

To solve for f, define the function H [0,«>)-*[0,°o) by 

H -1 (u) = ue u . (36) 

Note that H -1 is strictly increasing so it has an inverse, denot- 
ed H, which is also strictly increasing. Equation (34) can be 
written as 


f'/f = -(3/2) (t+t 0 ) -1 H(Kt(t+t 0 )" 2 f 2 ) 
and integrating while using (29) gives 


f(t) 


exp -(3/2) 

't 

(r+t 0 ) _1 H(Kr [r+t 0 ] _2 f 2 (r) ) dr 

L 

0 J 


(37) 
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Equation (37) can be solved for f by iteration 1 . 


Iteration of (37) is done numerically, but in numerical work we 
cannot go to the limit as t-*». We must work with a finite inter- 
val [0,T] and it is necessary to know how large T must be in 
order for f(T) to adequately approximate f(<»). This can be deter- 
mined by writing 


f (T)/f (») 


exp [( 3/2) 

00 

(r+t 0 )" 1 H(Kr[r+t 0 ]- 2 f 2 (r) ) dr 

L 

T J 


If we require, for example, that the ratio be between 0.99 and 1, 
then we require 


(3/2) 


(% 

00 

(r+t 0 )- 1 H(K7-[r+t 0 ] _2 f 2 (T)) 

JT 


dr < 0.01 . 


It is easy to show that H(u)<u for positive u. Also, f(r)<l so 
the above condition will be satisfied if 


1. A convergence proof is not required because convergence is not 
used to prove any theorems. Convergence is needed only for the 
numerical computation of f and can be verified on a case-by-case 
basis during the computation. But it is useful to note that (37) 

is similiar to (18) in the sense that upper bounds produce lower 
bounds and vice versa. If the initial guess is f=l (an upper 
bound) , iteration will produce a sequence of alternating upper 
and lower bounds. The correct solution is bracketed between any 
pair of adjacent iterates and is known to within a given toler- 
ance when the difference between adjacent iterates is less than 
that tolerance. 
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(3/2) 


00 


| Kr(T+t 0 ) 3 dr < 0.01 

JT 


or 


(3/2)K(T+t 0 /2) (T+t 0 ) -2 < 0.01 . 


This will be satisfied if 


T = 150K . (38) 

In addition to evaluating T, the previous analysis also demon- 
strated that the integral in (37) must be finite, even in the 
limit as t-oo. Therefore, we conclude from (37) that f(«)>0. It is 
not difficult to show that 

P E1 “*■ f ( t ) p 0 as t - oo (39) 

so that a lower bound N E1 (t) for the linear density N(t) satis- 
fies 


N Ei(t) -* f(t)N 0 as t -» oo (40) 

and does not go to zero as t-+oo. The conclusion is that some EHPs 
never will recombine (via AR) no matter how long we wait. It is 
not difficult to show that P £1 is within 1% of the right side of 
(39) and, therefore, N £1 is within 1% of the right side of (40) , 
when t>T . But f(t) is within 1% of f(°°) when t>T. Therefore, N E1 
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is within 2% of its limiting value f(°°)N 0 when t>T. Nearly all AR 
that ever will occur will do so in a time T given by (38) . 


With f now regarded as known, a lower bound estimate N £1 for 
the linear density can be evaluated from 


N E i(t) - 


p Ei d 2 X = 27 r 


00 

[2Ct+(fP 0 )“ 2 ] _1 / 2 r dr . (41) 


An upper bound for P is P £2 obtained from (18) according to 


P E2 = P CT C ( 47rD > 


-1 


P E i 3g d 2 X ' dt' 


(42) 


and an upper bound for the linear density is N E2 given by 


N E2( t ) - 


P E2 d X ‘ 


(43) 


Integrating (42) while using (21) and (28) gives 


N E2 (t) = N 0 -27rC 


/ 

{ 2Ct ' + [ f ( t ' ) p 0 (r,t')]~ 2 } ” 3 / 2 rdrdt'. (44) 


The integrals in (41) and (44) can be evaluated using (27) and 
the change in variables 


u = exp{ [4D(t+t 0 ) ] -1 r 2 ) 
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with the result 


N E1 

N E2 


fi" 1 ln[ f 2 + (l+f 2 2 ) 1/2 ] 


N 0 -(l/2) 


f 3 (f) df 


(45) 

(46) 


where 


f-L ■ (2Ct) 1/2 [47TD(t+t 0 ) ) _1 (47) 

f 2 = N^f (48) 

f 3 s (tf 1 ) _1 {ln[f 2 +(l+f 2 2 ) 1/2 ]-f 2 (l + f 2 2 )“ 1/2 ) . (49) 

4 . COMPUTER CODE AND NUMERICAL EXAMPLES 

A computer code that numerically evaluates the lower and upper 
bound linear density estimates N E1 and N E2 is listed in Appendix 
C. The value used for D is constructed from PISCES default values 
of 500 and 1000 cm 2 /V-s for hole and electron mobilities, respec- 
tively. This gives (at room temperature) D=17.3 cm /s. C is 
constructed from PISCES default values as discussed in section 2. 
The code assumes the initial density to be given by 


Pj = N 0 (7rr c 2 ) 1 exp(-r 2 /r c 2 ) 


(50) 


17 


where the initial linear density N 0 and characteristic radius r c 
are specified by the user. The code output is in a file called 
RCOM.OUT and is the estimated fractional surviving linear density 

[(N E i+N E2 )/ 2 ± (N E2 -N E1 )/2]/N 0 . 

Martin et al.[4] give the radial profile for a track produced 
by 270-MeV (initial energy) Kr at selected depths in a Si target. 
Their data was used to plot Figures 1, 2, and 3 which refer to 
depths of 5, 15, and 25 jim respectively. Because the computer 
code assumes a profile of the form (50) , it is necessary to curve 
fit. This motivates the construction of the upper curves in each 
figure, which are the integrated densities, i.e., the total 
number of EHPs contained within a given radius (the asymptotic 
values of these curves are the linear densities) . These curves 
are useful for telling us where curve fitting is and is not 
required to be accurate. Accuracy is not needed on intervals 
containing a negligible fraction of the EHPs. The dashed curves 
in the figures are the functions of the form (50) that were 
selected to represent the actual curves. In all cases, the fit- 
ting curves produce the same linear densities as the actual 
curves. The remaining criteria for selecting the fitting curves 
is a mixture of intuitive judgment and reasoning. For example, in 
the 25 urn case (Figure 3), the fit is reasonably good from 
r=5xl0 ^ to 5xl0 - ^ jim. This interval contains about two thirds 
of the total EHPs, and the fit is best where the density is 
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LOG-j Q(r/jim) 


Figure 1. Radial distribution of initial track density: Produced 
by 270-MeV (incident energy) Kr after traveling 5 /jm in Si. The 
lower solid curve is the volume density and the dashed curve is a 
fit to it. Upper solid curve is the total number of EHPs con- 
tained within the specified radius. 



LOG 10< P I cm J Q(10 10 cm-1) 
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Figure 2. Radial distribution of initial track density: Same as 
Figure 1 except that the depth is 15 /xm. 


0 


2 





Figure 3. Radial distribution of initial track density: Same 
Figure 1 except that the depth is 25 /im. 
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largest (with the exception of r<5xl0 4 /xm, but less than 10% of 
the EHPs are located there) which is also where AR is most impor- 
tant. 

The fitting curves in Figures 1, 2, and 3 are characterized by 
r c =1.5xl0 2 , l.lxlO -2 , and 2.76xl0 -3 /xm respectively, and 
Nq= 2 . 34xl0 10 , 2.58X10 10 , and 7.2xl0 9 cm -1 respectively. Using 
these numbers, the computer code predicts the fractional losses 
to be 5.0%, 10.4%, and 12.3% respectively. 


5. QUALITATIVE AND SEMI-QUANTITATIVE PREDICTIONS 

It was shown in section 3 that the linear density approaches a 
non-zero limiting value as t->a>. This behavior is due to diffusion 
which lowers the EHP density and, therefore, the AR rate (if D 
was set equal to zero in (1), we would find that the linear 
density goes to zero) . The AR rate decreases so fast with in- 
creasing t that some EHPs never will recombine (via AR) no matter 
how long we wait. It was also shown that the linear density is 
within approximately 2% of its asymptotic value after time T 
given by (38) (the word "approximately" is needed because it was 
the estimate P E1 rather than the actual P that was shown to be 
within 2% of its asymptotic value after time T) . Using numbers 
from the examples of section 4, we find that T is on the order of 
a ps . This is very fast compared to most charge transport pro- 
cesses, including the recovery stage of funnel ing [with the possi- 
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ble exception of funneling in a very thin (less than 5 micron) 
epi-layer] [5]. AR is not fast compared to the collapse stage of 
funneling [5], but the section of track below the collapsed 
depletion region remains intact (it is not pulled apart by the 
electric field) immediately after the collapse. If the track is 
dense enough to cause a nearly complete collapse (so that nearly 
all of the track remains intact, which is the expected situation 
for tracks dense enough for AR to be important) , the collapse 
stage will not significantly affect AR. Therefore, AR is expected 
to occur as section 3 predicted during a typical funneling pro- 
cess, and be essentially finished before any significant charge 
collection occurs. Therefore, when accompanying many charge 
transport processes, including funneling, AR should have the same 
effect as a reduction in ion linear energy transfer (LET) . We can 
define another LET, which will be called the "reduced LET", which 
relates to surviving linear density in the same way that ordinary 
LET relates to initial linear density. Charge collection calcula- 
tions should automatically include AR losses if reduced LET 
replaces ordinary LET. Single-event upset rate calculations 
should automatically include AR losses if experimentally measured 
device cross sections and Heinrich flux are plotted against 
reduced LET instead of ordinary LET. 

If we agree to call a 5 to 10% loss of EHPs "marginally impor- 
tant", then the examples of section 4 predict that AR is of 
marginal importance for 270-MeV Kr and, therefore, of major 
importance for ions exhibiting a significantly larger loss. 
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6. SOURCES OF ERROR 


There are a number of sources of error in these predictions of 
AR losses. The first is the replacement of a complicated system 
of nonlinear equations (the equations used by PISCES) with the 
simple diffusion equation (2) . The second source of error is the 
fact that the computer code (Appendix C) produces an approximate, 
not exact, solution to (2). The PISCES code is not practical for 
parametric studies (because of large CPU time) but it can be run 
once to access the errors stated above. A cylindrical coordinate 
version of PISCES was run for the 15 /zm example of section 4 . 
PISCES predicted a 10.1% loss compared to our prediction of a 
10.4% loss. Thus, it appears that the first two sources of error 
are not too serious. 

A third source of error is the need to curve fit a function of 
the form (50) to the initial EHP density, which might look noth- 
ing like a Gaussian. The amount of error depends on the individu- 
al case, but can be estimated by trying several fits and running 
the code for each to determine how sensitive the predicted AR 
loss is to Nq and r c . The error might be reduced by an analysis 
(which remains to be worked out) that tells us how to select good 
values for Nq and r c . 

The final source of error might be the most difficult to deal 
with. The equation used to solve for P assumes the charge carri- 
ers to be thermal ized (an assumption also made by PISCES) . But 
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significant AR requires r c to be on the order of 10 3 or 10 2 nm 
and, under these conditions, some of the AR occurs in a tiny 
fraction of a picosecond. It is not clear how thermalized the 
track is during AR. The impact of this source of error is unknown 
at this time and quantitative estimates should be regarded as 
unreliable until they are supported by experimental measurements. 


7. CONCLUSION 


Quantitative estimates should be regarded as unreliable at this 
time, as stated in section 6, but it is reasonable to expect the 
qualitative statements, regarding the existence and use of re- 
duced LET (see section 5) , to be valid. Although unreliable, the 
quantitative estimates give credibility to the assertion that AR 
is of marginal importance to 270-MeV Kr tracks and of major 
importance to denser tracks. Therefore, it would be highly de- 
sirable to design and conduct an experiment that can directly 
measure reduced LET for Kr and heavier ions. 
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APPENDIX A 


The objective is to show that iteration of (17) produces a 
sequence that converges to the solution P. It suffices to show 
that (17) satisfies the hypothesis of Banach's Fixed Point Theo- 
rem [6] . This will imply not only convergence but also that each 
iterate is at least as good of an approximation to P as the 
previous iterate . 

Write (17) as 


P = TP 


(Al) 


where T is the operator defined by 


TP = P B + (4jtD) 


-1 


[BP - CF (P) ] d 2 X' dt'. (A2) 


We will work with the set of continuous functions and the maximum 
norm, which produces a complete metric space [7]. T is obviously 
a self mapping, i.e., it maps continuous functions into continu- 
ous functions. The only remaining condition that must be estab- 
lished is the existence of an a<l such that 


| | TP 1 ~TP 2 | | < at | |P 1 -P 2 I I 


(A3) 


for all P 2 and P 2 . I 1*1 I denotes the maximum norm, i.e.. 
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| | f | | = max | f (r,t) | 
re [0,a>) 
te[0,T] 

where T is the largest value of t that the solution is required 
to apply to. To establish (A3) , note that 

[V(P 1 ) -V(P 2 ) ]e _B ^ t_t ' >G d 2 X' dt'll (A4) 


TP 1 -TP 2 I | = (47rD) 


-1 


where 


V(r) = Br-CF(r) = 3 CP m z t-CF(t) . 

Note that V is continuous and its derivative is defined and 
continuous except at two isolated points (r=±P M ). Therefore, 


V(P 1 )-V(P 2 ) = 


her 


(t) dr 


(A5) 


Except at the two isolated points, V' is bounded by 


V'(r) < 3CP m 2 . 


Combining this with (A5) and (A4) gives 
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||TP 1 -TP 2 || < (4ttD) 


-ll 


3CP m 2 |p 1 -P 2 | e ^ ^ G d 2 X 1 dt 1 


< 3CP m 2 (4ttD) 1 | |P 1 “P 2 




G d 2 X ' dt' II • (A6 ) 


The integral can be evaluated with the result 


-B(t-t') 


G d 2 X ' dt* = 47rD(3CP M 2 )" 1 [l-e Bt ] 


so that 


ft 


-B(t-t ' ) 


G d 2 X ' dt 1 


= max |4 ttD( 3CP m 2 ) 1 [l-e Bt ]| 
te[0,T] 

= 4ttD(3CP m 2 ) -1 [l-e _Bt ] 


and (A6) becomes 


where 


| TP-L-TP 2 | | < a | |P 1 “P 2 I I 


a = 1-e 


-BT 


Since a<l, the proof is complete. 
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APPENDIX B 


The objective is to prove (24) and (25). These statements are 
special cases of the more general statement which follows. Let 
j9 1 (X,t) and 0 2 (X,t) be bounded and well behaved (in the sense 
that div grad (3 and 6(3/6t exist) and let them satisfy 

/3 1 (X,0) = /3 2 (X,0) = P X (X) . (Bl) 


Define g ^ and g 2 by 

= C(3 ^ 3 - D div grad (3 ^ - SP^/St for i = 1,2 . (B2) 

The conclusion is that if g 2 (X, t) >g 1 (X, t) for all X and for all 
t>0, then p 2 (X,t)>p 1 (X,t) for all X and for all t>0. 


To prove this statement, let p B be a constant that is an upper 
bound for both \P-±\ and | P 2 1 • Define F by 


F*(T) 


T 3 if I T | < P B 

(3 b 2 t if \t\ > p B 


( B3 ) 


and write (B2) as 

D div grad p L - Sp L /St = CF *(P±) - g L (B4) 


which can be inverted (see section 3) to give 
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Pi = p B +( 47rD ) 


-1 


[B)0 i -CF*(^ i )+g i ]e -B ( t_t, )G d 2 X' dt • (B5) 


where 


B s 3C/3 b 2 


(B6) 


and P Q satisfies (12) and (13), except that B is now given by 
(B6) instead of (7) . The presence of g^ in (B5) does not affect 
the operators contraction mapping property and we conclude, from 
Appendix A, that (B5) can be solved by iteration, starting with 
any initial guess. For i=l,2, let (parentheses emphasize 
that j is a superscript, not an exponent) be defined by 


P 


i 


( 0 ) 


= 0 


( B7a) 


^(j + 1) = P B + 


(47TD) 


-1 


[B/3 i ^)-CF*(/3 i (3) )+ gi ]e -B ( t_t, >G d 2 X ' dt' 


( B7b) 


Note that 


Pi = lim Pi . 
j-+00 


Therefore, to show that P 2 -P p it suffices to show that 
p 2 ^^Pi^^ for all j =1 , 2 , . . . . This can be shown by induction. 
From 
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and (B7 ) we have p 2 ^ ^ 1 • Now assume that 




( B9 ) 


It suffices to show that this implies 


p 2 ^ +1 ) > 


(BIO) 


Note that Br-CF* ( r ) is an increasing function of r and assumption 
(B9) implies 


B/3 2 ^ ) -CF*(/3 2 (^) > B)3 1 ^)-CF*(/3 1 ^h . 


Combining this result with (B8) and (B7) produces (BIO) and 
completes the proof. 



APPENDIX C 


PROGRAM RECOM 
PARAMETER ( NT=2 0 1 ) 

DIMENSION T ( 0 : NT ) , F (NT) , F3 (NT) ,S(NT) 

REAL K,N,N1(NT) ,N2 (NT) 

COMMON T, F , K, TO 
C SET CONSTANTS . 

WRITE (*,*) 'ENTER INITIAL LINEAR TRACK DENSITY (1/CM) ' 

READ* , N 

WRITE (*,*) 'ENTER INITIAL CHARACTERISTIC RADIUS (CM)' 

READ* , RC 
D=17. 3 
C=3 . 79E— 31 
PI4=12. 56637 
PI4D=PI4 *D 
T0=RC*RC/(4.0*D) 

K=0 . 7357589*C*N*N/ (PI4D*PI4D) 

C CALL A SUBROUTINE TO CONSTRUCT THE TIME POINTS T ( 0) , . . , T (NT) . 

C ALL DATA TRANSFER IS THROUGH THE COMMON BLOCK. 

CALL TIME 

C CALL A SUBROUTINE TO CONSTRUCT THE ARRAY F WHICH REPRESENTS 

C THE FUNCTION f BY THE RELATION F ( J) =f (T ( J) ) . ALL DATA TRANSFER 

C IS THROUGH THE COMMON BLOCK. 

CALL FUN 

C CONSTRUCT AN ARRAY TO REPRESENT N1(T) GIVEN BY EQ.(43) AND AN 
C ARRAY TO REPRESENT F3 (T) GIVEN BY EQ.(47). 

DO 10 J=1 , NT 
TJ=T ( J) 

F1=SQRT ( 2 . 0*C*TJ) / (PI4D* (TJ+TO) ) 

F2=N*F1*F ( J) 

W1=SQRT ( 1 . 0+F2 *F2 ) 

W2=ALOG ( F2+W1) 

N1 ( J) =W2/F1 

F3 ( J) = (W2-F2/W1) / (TJ*F1) 

10 CONTINUE 

C CONSTRUCT THE INTEGRAL IN (44). S(J) IS THE INTEGRAL EVALUATED 
C AT t=T ( J) . CONSTRUCT N2(J) FROM THE INTEGRAL. 

S(1)=F3 (1) *T ( 1 ) 

N2 ( 1 ) =N-0 . 5*S ( 1 ) 

DO 20 J=1 , NT-1 

S (J+l) =0 .5*(F3(J)+F3(J+1) ) *(T(J+1)-T(J) )+S(J) 

N2 (J+l) =N-0. 5*S (J+l) 

20 CONTINUE 
C OUTPUT THE DATA. 

OPEN (UNIT=1 , FILE= ' RCOM. OUT ' , STATUS= ' UNKNOWN ' ) 

REWIND ( 1 ) 

WRITE ( 1 , * ) ' N = ' , N 
WRITE ( 1, *) ' RC= ' , RC 
WRITE (1,*)' ' 

WRITE ( 1, *) ' ' 

WRITE ( 1 , *) ' TIME (SEC) FRACTIONAL SURVIVING E-H PAIRS' 

DO 30 J=1 , NT 

AVG=0 . 5* (N1 ( J) +N2 (J) )/N 
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ERROR=0.5*(N2(J)-N1(J) )/N 
WRITE (1,40) T(J),AVG,'+ or ERROR 
30 CONTINUE 
CLOSE (1) 

40 FORMAT (El 2 . 4 , F16 . 3 , A8 , F6 . 3 ) 

END 

Q*** *************************************************************** * 

SUBROUTINE TIME 
PARAMETER ( NT= 201) 

DIMENSION T (0 : NT) , F (NT) 

REAL K 

COMMON T, F, K,T0 

C SELECT POINTS ON THE TIME AXIS AND DENOTE THEM AS T ( 0) , . . . T (NT) , 
C WHERE T (0) =0 , T(1)=T0/100, T(NT)=150K, AND THE REMAINING Ts 
C ARE CHOSEN SO THAT T ( 1) , . . . , T (NT) ARE EQUALLY SPACED ON A 
C LOGARITHMIC SCALE. 

T (0) =0 . 0 
T ( 1) =T0/100 . 0 
T (NT) =150 . 0*K 
RNT=FLOAT (NT) 

A=ALOG10 (T (NT) /T ( 1) )/(RNT-1.0) 

B=ALOG10 (T ( 1) /TO) -A 
DO 10 J=2 , NT-1 
RJ=FLOAT ( J) 

T ( J) =T0*10 .0** (A*RJ+B) 

10 CONTINUE 
RETURN 
END 

Q** **************************************************************** * 

SUBROUTINE FUN 
PARAMETER ( NT= 201) 

DIMENSION T (0 : NT) , F (NT) ,H1(NT) ,S(NT) ,F1(NT) 

REAL K 

COMMON T, F, K, TO 

C THE FUNCTION f IN EQ.(37) IS REPRESENTED BY THE ARRAY F AND IS 
C SOLVED BY ITERATION. FIRST CONSTRUCT THE INITIAL GUESS. 

DO 10 J=1 , NT 
F ( J) =1 . 0 
10 CONTINUE 

C START THE ITERATION LOOP. 

20 CONTINUE 

C CONSTRUCT ARRAYS TO SIMPLIFY THE INTEGRAL IN (37) . HI IS THE 
C INTEGRAND AND W IS THE ARGUMENT OF H. H IS A FUNCTION 
C SUBPROGRAM. 

DO 30 J=1 , NT 

W=K*T ( J) *F ( J) *F(J)/ (T(J) +T0) **2.0 
HI ( J) =1 . 5*H (W)/(T(J) +T0) 

30 CONTINUE 

C CONSTRUCT THE INTEGRAL IN (37). S(J) IS THE INTEGRAL 
C EVALUATED AT T(J). 

S ( 1) =0 . 5*H1 ( 1) *T ( 1) 

DO 40 J=1 , NT-1 

S ( J+l) =S (J)+0.5*(H1(J) +H1 ( J+l) ) *(T(J+1) -T(J) ) 

40 CONTINUE 

C STORE THE EXISTING Fs IN A NEW ARRAY, FI, FOR A FUTURE 
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C ACCURACY CHECK, AND UPDATE THE VALUES STORED IN F USING (37) . 

DO 50 J=1 , NT 
F1(J)=F(J) 

F(J)=EXP(-1.0*S(J) ) 

50 CONTINUE 

C PERFORM ACCURACY CHECK AND ESCAPE FROM LOOP IF SATISFIED. 
IERROR=0 
DO 60 J=1 , NT 

IF (ABS(F1(J) -F(J) ) .GE. 0.0001) IERROR=l 
60 CONTINUE 

IF (IERROR.EQ. 1) GO TO 20 
DO 70 J=1 , NT 
F(J)=0.5*(F(J) +F1 ( J) ) 

70 CONTINUE 
RETURN 
END 

C************************************************************* 
FUNCTION H (X) 

C THE OBJECTIVE IS TO SOLVE FOR U=H(X). THIS IS DONE BY ITERATING 
C THE EQUATION U=0 . 5* (LN (X/U) +U) IF X >= .9 OR THE EQUATION 
C U=X*EXP(-U) IF X < .9. USE THE INITIAL GUESS U=X/2. 

U=0 . 5*X 

10 IF (X.GE.0.9) THEN 

U=0 . 5* (ALOG (X/U) +U) 

ELSE 

U=X*EXP(-U) 

END IF 

IF (ABS (U*EXP(U)/X-1. 0) .LE. 0. 0001) THEN 
H=U 

GO TO 20 
END IF 
GO TO 10 
20 CONTINUE 
RETURN 
END 
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